DALS021-MDS与PCA

发表于 | 分类于
| 字数统计: 2,527 | 阅读时长 ≈ 11

前言

这一部分是《Data Analysis for the life sciences》的第8章统计模型的第2小节，这一部分的主要内容涉及MDS和PCA，相应的Rmarkdown文档可以参考作者的Github。

MDS

MDS的全称为multi-dimensional scaling，即多维数据缩放。在这 一部分中，我们会使用基因表达的数据来作为案例讲解一下。为了简化说明，我们仅考虑3个组织：

1
2
3
4
5
6
7
library(rafalib)
library(tissuesGeneExpression)
data(tissuesGeneExpression)
colind <- tissue%in%c("kidney","colon","liver")
mat <- e[,colind]
group <- factor(tissue[colind])
dim(mat)

结果如下所示：

1
2
> dim(mat)
[1] 22215    99

现在我们要研究一下这个数据集，我们想知道，储存在mat列中的基因表达谱的数据在不同的组织间的相似性如何。由于数据很大，无法直接画出相应的多维点图。我们通常只能绘制出二维图形，如果我们要绘制出每两个样本之间的基因表达情况不现实。而MDS图形就是为了解决这个问题而提出来的。

MDS背后的数学原理

前面我们已经知道了SVD和矩阵代数，那么我们理解MDS就相对清楚了。为了说明MDS，我们先来看一下SVD分解，如下所示：

我们假设 $\mathbf{U^\top Y=DV^\top}$ 的前两列的平方和剩余列的平方和。因此它们可以写为$d_1+ d_2 \gg d_3 + \dots + d_n$ 其中 $d_i$ 是$\mathbf{D}$ 是第i列（原文是i-th entry）。当出现这种情况时，我们就会得到如下公式：

这就表明，第$i$列近似等于：

如果我们们定义下面的二维向量：

那么：

上面的这个推导告诉我们，在样本$i$和样本$j$之最的距离近拟等于下面二维数据点的距离：

因为$Z$是一个二维向量，因此我们可以通过绘制$\mathbf{Z_{1}}$和$\mathbf{Z_{2}}$来发展示这两个样本的距离。现在我们绘制出它们的距离：

1
2
3
4
5
6
s <- svd(mat-rowMeans(mat))
PC1 <- s$d[1]*s$v[,1]
PC2 <- s$d[2]*s$v[,2]
mypar(1,1)
plot(PC1,PC2,pch=21,bg=as.numeric(group))
legend("bottomright",levels(group),col=seq(along=levels(group)),pch=15,cex=1.5)

从图片上我们可以看出，数据点按照相应的组织区分开来了。上面的这种分开的精确近似取决于前两个主成分解释变异的程度。像上面那样所示，我们可以绘制出每个主成分可以解释的变异程度：

1
plot(s$d^2/sum(s$d^2))

虽然前两个主成分解释了超过50%的变异，不过前面的图形还是没有展示出大量的信息。但是这种图已经足够用于进行可视化大量的数据了。此外，我们还可以注意到，我们能够绘制其它的主成分来研究这些数据点，例如我们绘制第3个和第4个主成分：

1
2
3
4
5
PC3 <- s$d[3]*s$v[,3]
PC4 <- s$d[4]*s$v[,4]
mypar(1,1)
plot(PC3,PC4,pch=21,bg=as.numeric(group))
legend("bottomright",levels(group),col=seq(along=levels(group)),pch=15,cex=1.5)

从上面图形中我们可以看到，第4个主成分能够将肾脏组织的样本强烈分开。在后面的部分中，我们会讲到批次效应(batch effects)会解释这种情况。

cmdscale()函数

我们在上面使用了svd()函数来进行计算，不过R中有一个专门的函数用于计算MDS，生成MDS图。这个函数就是cmdscale()函数，这个函数将距离对象作为参数，然后使用主成分分析来对这些距离进行近似计算。这个函数比使用svd()函数更高效（因为不可能实现完全的svd()函数计算，那样比较花时间）。此函数默认返回二维的数据，不过我们通过设定参数k（默认情况下，k=2）可以改变结果中的维度：

1
2
3
4
5
6
d <- dist(t(mat))
mds <- cmdscale(d)
mypar()
plot(mds[,1],mds[,2],bg=as.numeric(group),pch=21,
xlab="First dimension",ylab="Second dimension")
legend("bottomleft",levels(group),col=seq(along=levels(group)),pch=15)

再看另外一个：

1
2
3
4
5
6
mypar(1,2)
for(i in 1:2){
plot(mds[,i],s$d[i]*s$v[,i],main=paste("PC",i))
b = ifelse( cor(mds[,i],s$v[,i]) > 0, 1, -1)
abline(0,b) ##b is 1 or -1 depending on the arbitrary sign "flip"
}

任意符号

SVD并非是唯一的，只要我们用-1乘以$\mathbf{U}$的样本列，我们就能使用-1乘以$\mathbf{V}$的任意列，通过下面的转换我们就能看出来（这一段不懂）：

扣除平均值

在所有的计算中，当我们计算SVD时，都会扣除行(row)的均值。如果我们要试图计算两列之间的近似距离，那么在$\mathbf{Y}_{i}$和$\mathbf{Y}_{j}$之间的距离就与$\mathbf{Y}_i - \mathbf{\mu}$和$\mathbf{Y}_j - \mathbf{\mu}$之间的距离相同，因为当我们过计算时，中间的$\mu$就会被消去：

因为扣除行均值可以降低总的变异，它可以使得SVD的结果近更为逼近。

练习

P357

PCA

PCA的相关资料可以参考作者的Github。

前面我们已经提到了PCA，这里继续深入一步，讲一下PCA背后的数学原理。

案例：双胞胎身高

我们先使用模拟数据的案例展示一个旋转，这个旋转与PCA有着很大的有关系：

1
2
3
4
5
6
7
8
9
10
library(rafalib)
library(MASS)
n <- 100
set.seed(1)
Y=t(mvrnorm(n,c(0,0), matrix(c(1,0.95,0.95,1),2,2)))
mypar()
thelim <- c(-3,3)
plot(Y[1,], Y[2,], xlab="Twin 1 (standardized height)", 
     ylab="Twin 2 (standardized height)", xlim=thelim, ylim=thelim)
points(Y[1,1:2], Y[2,1:2], col=2, pch=16)

这里我们专门来解释一下什么是什么成分（principla components）。

我们使用 $\mathbf{Y}$ 这个 $2 \times N$ 矩阵来表示我们的数据。这个类似于我们检测了两组基因的信息，每列表示1个样本。现在我们的任何就是，找到一个 $2 \times 1$ 向量 $\mathbf{u}_1$ ，使其满足 $\mathbf{u}_1^\top \mathbf{v}_1 = 1$，它能使 $(\mathbf{u}_1^\top\mathbf{Y})^\top (\mathbf{u}_1^\top\mathbf{Y})$ 最大。这个过程可以被视为每个样本，或$\mathbf{Y}$向子空间 $\mathbf{u}_1$ 的投影。因此，我们需要将坐标系进行置换，使新的坐标系能够显示出最大变异。

我先试一下 $\mathbf{u}=(1,0)^\top$。这个投影公仅能够给出双胞胎1的身高（橘黄色）。图片标题中显示的是平方和。

projection_not_PC1, fig.align
1
2
3
4
5
6
mypar(1,1)
plot(t(Y), xlim=thelim, ylim=thelim,
     main=paste("Sum of squares :",round(crossprod(Y[1,]),1)))
abline(h=0)
apply(Y,2,function(y) segments(y[1],0,y[1],y[2],lty=2))
points(Y[1,],rep(0,ncol(Y)),col=2,pch=16,cex=0.75)

我们能否找到一个方向，使得坐标系旋转后，能够表示更高的变异？例如

$\mathbf{u} =\begin{pmatrix}1\-1\end{pmatrix}$ 这个怎么样？它不满足 $\mathbf{u}^\top\mathbf{u}= 1$ ，因此我们可以使用另外一个向量，即
$\mathbf{u} =\begin{pmatrix}1/\sqrt{2}\-1/\sqrt{2}\end{pmatrix}$

projection_not_PC1_either, fig.cap
1
2
3
4
5
6
7
8
9
10
11
12
13
library(rafalib)
u <- matrix(c(1,-1)/sqrt(2),ncol=1)
w=t(u)%*%Y
mypar(1,1)
plot(t(Y),
     main=paste("Sum of squares:",round(tcrossprod(w),1)),xlim=thelim,ylim=thelim)
abline(h=0,lty=2)
abline(v=0,lty=2)
abline(0,-1,col=2)
Z = u%*%w
for(i in seq(along=w))
  segments(Z[1,i],Z[2,i],Y[1,i],Y[2,i],lty=2)
points(t(Z), col=2, pch=16, cex=0.5)

这个图形与双胞胎的差异有关，我们知道这个差异很少的。通常平方和我们可以确实这一点，最后我们试一下这个向量：

PC1, fig.cap
1
2
3
4
5
6
7
8
9
10
11
12
13
u <- matrix(c(1,1)/sqrt(2),ncol=1)
w=t(u)%*%Y
mypar()
plot(t(Y), main=paste("Sum of squares:",round(tcrossprod(w),1)),
     xlim=thelim, ylim=thelim)
abline(h=0,lty=2)
abline(v=0,lty=2)
abline(0,1,col=2)
points(u%*%w, col=2, pch=16, cex=1)
Z = u%*%w
for(i in seq(along=w))
  segments(Z[1,i], Z[2,i], Y[1,i], Y[2,i], lty=2)
points(t(Z),col=2,pch=16,cex=0.5)

这个图形与重新缩放(re-scaled)后的平均高度有关，它有着最大的平方和。这是一个数学计算程序，它能够计算出一个 $\mathbf{v}$ ，能够使平方和最大，SVD就是这样的一个程序。

主成分

正交向量能够使平方和最大：

$\mathbf{u}_1^\top\mathbf{Y}$ 指的就是第1PC。e用于获得PC的加权(weights) $\mathbf{u}$ 指的就是因子载荷(loadings)。使用旋转这种操作，它指的就是第1PC的旋转方向。

为了获得第2PC，我们可以重复上述操作，但是残差如下：

第2PC的向量含有以下性质：

它能使 $(\mathbf{rv}_2)^\top \mathbf{rv}_2$最大，

当 $Y$ 是 $N \times m$ 时，我们可以重复地找到第3，第4，第5，等主成分。

prcomp

我们已经介绍了如何使用SVD来计算PC。介理，R中有一个专门的函数可以用于找到主成分，即prcomp()，在这个案例中，数据默认中心化的，这个函数的使用如下所示：

1
pc <- prcomp( t(Y) )

计算出的结果与SVD相同，直到符号翻转（produces the same results as the SVD up to arbitrary sign flips，实在没理解这句话什么意思）

pca_svd, fig.cap
1
2
3
4
5
s <- svd( Y - rowMeans(Y) )
mypar(1,2)
for(i in 1:nrow(Y) ){
  plot(pc$x[,i], s$d[i]*s$v[,i])
}

因子载荷可以通过下面方式计算：

1
pc$rotation

计算结果如下所示：

1
2
3
4
> pc$rotation
           PC1        PC2
[1,] 0.7072304  0.7069831
[2,] 0.7069831 -0.7072304

它就相当于 (up to a sign flip？这个不懂) :

1
s$u

计算结果如下所示：

1
2
3
4
> s$u
           [,1]       [,2]
[1,] -0.7072304 -0.7069831
[2,] -0.7069831  0.7072304

解释的方差等价于：

1
pc$sdev

计算结果如下所示：

1
2
> pc$sdev
[1] 1.2542672 0.2141882

现在我们将Y转置一下，因为prcomp()函数与我们平时所用的高通量数据储存有点不太一样，平时我们的数据是列为样本，行为特征值，而prcomp()函数则是正好相反，它处理的数据列是特征值，行是样本名。